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Abstract 

In [N. Ujevic, New iterative method for solving linear systems, Appl. Math. Comput. 
179 (2006) 725730], a new iterative method for solving linear system of equations was 
presented which can be considered as a modification of the Gauss-Seidel method. 
Then in [Y.-F. Jing and T.-Z. Huang, On a new iterative method for solving linear 
systems and comparison results, J. Comput. Appl. Math., In press] a different 
approach, say 2D-DSPM, and more effective one was introduced. In this paper, we 
improve this method and give a generalization of it. Convergence properties of this 
kind of generalization are also discussed. We finally give some numerical experiments 
to show the efficiency of the method and compare with 2D-DSPM. 
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1. Introduction 

Consider the linear system of equations 

Ax = b, (1) 

where A e M. nxn is a symmetric positive definite (SPD) matrix and x, b G M. n . The 
Gauss-Seidel method is an stationary iterative method for solving linear system of 
equation and is convergent for SPD matrices. This method is frequently used in 
science and engineering, both for solving linear system of equations and precondi- 
tioning [21 EJ. It can be easily seen that the Gauss-Seidel method is an special case of 
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a projection method [HE]- Let /C and £ be two m-dimensional subspaces of M n . Let 
also x be an initial guess of the solution of ([T]). A projection method onto fC and 
orthogonal to £ is a process which finds an approximate solution i G K n to |I| by- 
imposing the conditions that x e io + /C and the new residual vector be orthogonal 
to C (Petrov-Galerkin condition), i.e. 

Find x G xq + /C, such that b — AxJ-C. (2) 

It is well known that an iteration of the elementary Gauss-Seidel method can be 
viewed as a set of projection methods with C = fC = {e^}, % = 1,2, ... ,n, where e, 
is the ith column of the identity matrix. In fact, a single correction is made at each 
step of these projection steps cycled for i — 1, . . . , n. 

In [7], Ujevic proposed a modification of the Gauss-Seidel method which may 
be named as a "one-dimensional double successive projection method" and referred 
to as 1D-DSPM. In an iteration of 1D-DSPM, a set of double successive projection 
methods with two pairs of one-dimensional subspaces are used. In fact, in an iter- 
ation of 1D-DSPM two pairs of subspaces ()Ci,£i) and (/C2, £2) of one dimension 
are chosen while it makes double correction at each step of the process cycled for 
% = 1, 2, . . . , n. In jl], Jing and Huang proposed the "two-dimensional double succes- 
sive projection method" and referred to as 2D-DSPM. In an iteration of 2D-DSPM, 
a set of projection methods with a pairs of two-dimensional subspaces /C and C is 
used and a double correction at each step of the projection steps is made. 

In this paper, a generalization of 2D-DSPM say mD-SPM which is referred to 
as "m-dimensional successive projection method" is proposed and its convergence 
properties are studied. For m = 2, the mD-SPM results in 2D-DSPM. 

Throughout this paper we use some notations as follows. By (., .) we denote the 
standard inner product in MJ 1 . In fact, for two vectors x and y in M. n , (x,y) = y T x. 
For any SPD matrix M G M nxn , the M-inner product is defined as (x, jj)m — (Mx, y) 
and its corresponding norm is ||o;||m = (x,x) 1 ^. 

This paper is organized as follows. In section 2, a brief description of 1D-DSPM 
and 2D-DSPM are given. In section 3, the mD-SPM is presented and its convergence 
properties are studied. In section 4 the new algorithm and its practical implemen- 
tations are given. Section 5 is devoted to some numerical experiments to show the 
efficiency of the method and comparing with 1D-DSPM and 2D-DSPM. Some con- 



2 



eluding remarks are given in 6. 



2. A brief description of 1D-DSPM and 2D-DSPM 

We review 1D-DSPM and 2D-DSPM in the literature of the projection methods. 
As we mentioned in the previous section in each iteration of 1D-DSPM a set of double 
successive projection method is used. Let Xk be the current approximate solution. 
Then the double successive projection method is applied as following. The first step 
is to choose two pairs of the subspaces /Ci = C\ — {v\}, IC 2 = C 2 = {v 2 } and the 
next approximate solution Xk+i is computed as follows 

Find Xk+i G Xfc + JCi, such that b — Axk+iJ-Ci, (3) 

Find Xk+i € Xk+i + JC 2 , such that b — Axk+i-i-C 2 . (4) 
This framework results in [Jj, [7j 

x fc+ i = x k + axvt + (3 2 v 2 

where 

"i = -Pi/a, P2 = (epi - ap 2 )/ad, (5) 

in which 

a=(v 1 ,v 1 ) A , c=(v 1 ,v 2 ) A , d=(v 2 ,v 2 ) A - (6) 

In [7], it has been proven that this method is convergent to the exact solution of 
(JTj) for any initial guess. 

In the 2D-DSPM, two two-dimensional subspaces K = C = span{v!,v 2 } are 
chosen and a projection process onto /C and orthogonal to C is used instead of dou- 
ble successive projection method used in 1D-DSPM. In other words, two subspaces 
fC — C — span{t>i, t>2} are chosen and a projection method is defined as following. 

Find Xk+i € Xk + /C, such that b — Axk+i-i-C (7) 
In [I], it has been shown that this projection process gives 

a _ C P 2 ~ ^ _ epi - ap 2 

ad — c 2 ' ad — c 2 ' 
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where pi, P2, a, b, and c were defined in Eqs. fl5]) and (jH])- It has been proven in 
[I] that the 2D-DSPM is also convergent. Theoretical analysis and numerical exper- 
iments presented in [4] show that 2D-DSPM is more effective than the 1D-DSPM. 

A main problem with 1D-DSPM and 2D-DSPM is to choose the optimal vec- 
tors Vi and i>2- in this paper, we first propose a generalization of 2D-DSPM and then 
give a strategy to choose vectors v± and v 2 in a special case. 



3. m-dimensional successive projection method 

Let {ui,i>2, • • • ,v m } be a set of m independent vectors in W 1 . For later use, let 
also V m = [v 1, t>2, . . • , v m ]. Now we define the mD-SPM as follows. In an iteration of 
the mD-SPM we use a set of projection process onto K = span{t>!, t> 2 , . . . ,v m } and 
orthogonal to C = JC. In other words, two m-dimensional subspaces fC and C are used 
in the projection step instead of two two-dimensional subspaces used in 2D-DSPM. 
In this case Eq. (j2J) turns the form 

Find y m G 1R" 1 such that x k+l = x k + V m y m , and V^(b - Ax k+ i) = 0. (8) 



We have 



r k +i 



b - Ax k+1 

b - A(x k + V m y m ) 

r k - AV m y 

in ■ 



where r k = b — Ax k . Hence from Eq 



(jH) we deduce 



= V%r k+1 = VZ(r k - AV m y) 



Vlr k - VlAV mVl 



=> V^AV m y m = V£r k . 



The matrix V^AV m is an SPD matrix, since A is SPD. Therefore 



Um 



(k: 



ZAV m y l V^r h . 



(9) 



Hence, from (jSj) we conclude that 



x k+1 =x k + V m (V^AV m )- l V^r k . 



(10) 
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Theorem 1. Let A be an SPD matrix and assume that Xk is an approximate solution 
o/([ip. Then 

IMfclU ^ IMfc+ilU) (11) 

where dk = x* — Xk and dk+i = x* — Xk+i in which Xk+i is the approximate solution 
computed by Eq. [Tty) . 

Proof. It can be easily verified that dk+i = dk — V m y m and Adk = r k where y m is 
defined by ([7]). Then 

(Ad k+ i, d k+ i) = {Ad k - AV m y m} d k - V m y m ) 

= (Ad k , d k ) - 2(V m y m , r k ) + (AV m y m , V m y m ) by Ad k = r k 
= (Ad k , d k ) - (y m , Vjjk) by ®. 

Therefore 

IMfclll — = (Ad k ,dk) — (Adk+i, d k+ i) 

= (Vm, V*r k ) 

= (V^rkfiV^AVj-'V^n, by© 
=: S(r k ). 

Since (V^AVm)' 1 is SPD then we have 

IM&IIa ~ IMfc+ilIyi — 0) 
and the desired result is obtained. □ 

This theorem shows that if V^r k = then S(rk) = and we don't have any 
reduction in the square of the A-norm of error. But, if V^r k ^ then the square of 
the A-norm of error is reduced by S(r k ) > 0). 

Theorem 2. Assume that A is an SPD matrix and £ = JC. Then a vector Xk+i 
is the result of projection method onto K orthogonal to C with the starting vector Xk 
iff' it minimizes the A-norm of the error over Xk + /C. 

Proof. See g], page 126. □ 

This theorem shows that if /Ci = L\ C /C2 = £2, then the reduction of A-norm of 
the errors obtained by the subspaces /C2 = £2 is more than or equal to that of the 
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subspaces /Ci = C\. Hence by increasing the value of m, the convergence rate may 
increase. 



In the continue we consider the special case that the vectors Vi are the column 
vectors of the identity matrix. The next theorem not only proves the convergence of 
the method in this special case but also gives an idea to choose the optimal vectors v j. 

Theorem 3. Let 12, . . . , i m } be the set of indices of m components of largest 
absolute values in r k such that %\ < %i < . . . < i m . If Vj = e^, j = 1, . . . , m then 

Til 

\\d\\ A -\\d new \\ A > \\r k \\l (12) 

Proof. Let E m = [e^, e i2 , . . . , e im }. By Theorem 1, we have 



ll^fcll^ — IMfc+illi — S( r k) 

= {Elr k ) T {ElAE m )- l Elr k . 

Then by using Theorem 1.19 in [6J we conclude 

S(r k ) > ^{{ElAE^WElnWl > 1 A \E T m r k \\l (13) 

where for a square matrix Z, X min (Z) and A max (Z) stand for the smallest and largest 
eigenvalues of Z. It can be easily verified that [3] 

{r k ,r k ) n 



Hence 



S[Tk) £ OT)l |r ' £ (15) 



and the desired result is obtained. □ 

Eq. (12) shows the convergence of the method. Eq. (fT3"|) together with the 
first relation of the equation (TTJJ give 

S(r„) > —L— ||*£r fc ||2. 
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This equation gives an idea to choose indices ij, j = 1, . . . , m. In fact, it shows that 
if these indices are the m components of the largest absolute values in r^, then the 
lower bound of S{rk) depends on ||£'^rfc|||, and will be as large as possible. 

In [3], another theorem for the convergence of the method obtained by vi = was 
presented and an algorithm based upon this theorem was constructed for computing a 
sparse approximate inverse factor of an SPD matrix and was used as a preconditioner 
for SPD linear system. 

4. Algorithm and its practical implementations 

Hence, according to the results obtained in the previous section we can summa- 
rized the mD-DSM in the special case that following:. 

Algorithm 1: mD-DSM 

1. Choose an initial guess x G lR n to ([I]) and r = b — Ax . 

2. Until convergence, Do 



3. x := Xo 

4. For i — 1, . . . , n, Do 

5. Select the indices ii, i 2 , ■ ■ ■ , i m of r as defined in Theorem 3 

6 . E m ■ [cjj , e i 2 , . . . , e j m ] 

7. Solve (E^AE m )y m = E^r for y m 

8. x . x ~\~ E rn y m 

9. r := r - AE m y m 

10. EndDo 

11. Xq := x and if x Q has converged then Stop 



12. EndDo 

In practice, we see that the matrix E^AE m is a principal submatrix of A with 
column and row indices in {i±,i2, ■ ■ ■ ,i m }- Hence, we do not need any computation 
for computing the matrix E^AEm in step 7. For solving the linear system in step 
7 of this algorithm one can use the Cholesky factorization of the coefficient matrix. 
Step 8 of the algorithm may be written as 
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• For j = 1, 2, . . . , m 

• EndDo 

Hence only m components of the vector x are modified. Step 9 of the algorithm can 
be written as 

r = r- ^2(y m )j (kj 

3=1 

where is the ij column of the matrix A. 

It can be seen that Algorithm 2 in [1] is an special case of this algorithm. In 
fact, if m = 2 and the indices i\ and 12 are chosen as i\ = i and 12 = i — ij ga p 
{12 = i — ijgap + n if i < ijgap), where ij gap is a positive integer parameter less than 
n, then Algorithm 2 in [4J is obtained. 

As we see, the first advantage of our algorithm over Algorithm 2 in [1] is that 
our algorithm chooses the indices {11,12, ■ ■ ■ ,i m }, automatically. Another advantage 
is that our algorithm chooses the indices such that the reduction in the square of the 
A-norm of the error is more than that of Algorithm 2 in [4| . Numerical experiments 
in the next section also confirm also this fact. The main advantage of Algorithm 2 
in [1] over our algorithm is that only m components of the current residual are com- 
puted whereas in our algorithm the residual vector should be computed for choosing 
m indices of largest components in absolute value. 



5. Numerical experiments 

In this section we give some numerical experiments to compare our method with 
Algorithm 2 in [3]. Numerical results have been obtained by some MATLAB codes. 
We use all of the assumptions such as initial guess, exact solution, stopping criterion, 
and the examples used in [3]. Let b = Ae, where e is an n- vector whose elements 
are all equal to unity, i.e., e = (1, 1, . . . , 1) T . We use || a^^+i — ^fc||oo < 10~ 6 as the 
stopping criterion. An initial guess equal to ), where x l = 0.001 x i, 

z = l,2,...,nis chosen. For each of the systems we give the numerical experiments 
of Algorithm 2 in [1] with ij gap = 2 and 500 and our algorithm with m = 2, 3, 4 and 5. 
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lable 1: Results 


for the Example 1 




Algorithm 2 in [4J 2D-DSM 


3D-DSM 4D-DSM 


5D-DSM 


6 (ijgap = 2) 5 


4 3 


2 


7 feap = 500) 






Table 2: Results 


for the Example 2 




Algorithm 2 in [4J 2D-DSM 


3D-DSM 4D-DSM 


5D-DSM 


8 (tj gap = 2) 7 


6 4 


4 


9 fea P = 500) 







Example 1. Let A = (a^) where 

an = 4n, a M+ i = a i+M = n, =0.5 for i = 1, 2, . . . , n, j ^ i, i + 1. 

We also assume n = 1000. Numerical experiments in terms of iteration number were 
shown in Table 1. 

Numerical experiments presented in Table 1 show that the 2D-DSM method gives 
better results than the Algorithm 2 in [4J . This table also shows the effect of increase 
in m on the number of iterations for convergence. 

Example 2. Let A = (a^-) be the same matrix used in the previous example except 
the diagonal entries are changed to 

an = 3n, i — 1, 2, . . . , n. 

Numerical experiments were given in Table 2. 

This table also shows the advantages of our method on Algorithm 2 [4]. 

Example 3. Our third set of test matrices used arise from standard five point 
finite difference scheme to discretize 

-An + a(x,y)u x + b(x,y)u y + c(x,y)u = f(x,y), in f2 = [0, 1] x [0, 1], 

where a(x, y), b(x, y), c(x, y) and d(x, y) are given real valued functions. We consider 
three following cases: 
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Table 3: Results for the Example 3 



Cases Algorithm 2 in [4j 2D-DSM 3D-DSM 4D-DSM 5D-DSM 



Case 1 


391 
323 


{ijgap — 2) 
{ijgap = 500) 


226 


153 


116 


94 


Case 2 


312 
256 


{ijgap = 2) 
{ijgap = 500) 


192 


131 


100 


80 


Case 3 


302 
250 


{ijgap = 2) 
{ijgap = 500) 


218 


151 


115 


93 



Casel : a(x, y) = 0, b(x, y) = 10(x + y), c(x, y) = 10{x - y), f(x,y)=0, 
Case2 : a(x, y) = -10(x + y), b(x,y) = -10(x - y), c(x,y) = l, f(x,y) = 0, 
Case3 : a(x, y) = 10e xy , b(x,y) = I0e~ xy , c(x,y) — 0, f(x,y) — 0. 
We assume m = 32. In this case we obtain three SPD matrices of order n = 32 x 32 
[1] and used them as the coefficient of the linear systems. Numerical results were 
given in Table 3. 

This table also confirm that our method is more effective that the Algorithm 2 

i- 



5. Conclusion 

In this paper a generalization of the 2D-DSPM [1] which itself is a generalization 
of 1D-DSPM |7J is presented. 1D-DSPM and 2D-DPM need to prescribed some sub- 
spaces of lR n for the projection steps. But our method in the spacial case chooses this 
subspaces automatically. Theoretical analysis and numerical experiments presented 
in this paper showed that our method is more effective that 2D-DSPM. 
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